Wnioski końcowe

  1. Głównym czynnikiem wpływającym na długość śledzia jest temperatura przy powierzchni wody.
  2. Zagęszczenie planktonu jest skorelowane z długością śledzia w niewielkim stopniu. Jego spadek jest prawdobnie także następstwem wzrostu temperatury wody.
  3. Intensywność połowów ma niewielki wpływ na przeciętną długość śledzia.

Wprowadzenie

Niniejszy raport stanowi ćwiczenie umiejętności oraz prezentację wiadomości zdobytych podczas przedmiotu Programowanie w R i innych zajęć w ramach studium podyplomowego __Hurtowie Danych i analiza danych w celach biznesowych_ Politechniki Poznańskiej. Przy okazji podejmiemy próbę odpowiedzi na pytanie postawione w zadaniu. Co powoduje karłowancenie śledzi oceanicznych wyławianych w europie?

Dane, które posłużyły do badania można pobrać tutaj.

W celu odpowiedzi na postawione pytanie posłużymy się narzędziami exploracji danych zaimplementowanymi w języku R. Wykorzystane biblioteki prezentuje poniższy bloku kodu.

library(gridExtra) 
library(reshape2) 
library(dplyr)
library(ggplot2)
library(knitr)
library(corrplot) 
library(plotly)
library(caret)

Dane źródłowe

Surowe dane są wynikiem pomiarów rozmiaru śledzia oceanicznego wyławianego w Europie. Dane były zbierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi. Ponadto załączona do pliku informacja świadczy, że obserwacje były zapisywane chronologicznie. W poniższej tabeli zamieściłem opis atrybutów źródłowego zbioru danych. Tam gdzie było to możliwe załączyłem linki do odpowiednich artykułów wikipedii.

nazwa zmiennej opis [jednostka]
X Liczba porządkowa
length długość złowionego śledzia [cm]
cfin1 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]
cfin2 dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]
chel1 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]
chel2 dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]
lcop1 dostępność planktonu [zagęszczenie widłonogów gat. 1]
lcop2 dostępność planktonu [zagęszczenie widłonogów gat. 2]
fbar natężenie połowów w regionie [ułamek pozostawionego narybku]
recr roczny narybek [liczba śledzi]
cumf łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]
totaln łączna liczba ryb złowionych w ramach połowu [liczba śledzi]
sst temperatura przy powierzchni wody [°C]
sal poziom zasolenia wody [Knudsen ppt]
xmonth miesiąc połowu [numer miesiąca]
nao oscylacja północnoatlantycka [mb]

Opis zbioru przeznaczonego do analizy nie zawiera informacji na temat metodologii pomiaru zagęszczenia planktonu. Nie znamy też celu, dla którego zdecydowano się na pomiar dwóch szczególnych gatunków widłonogów, ani dlaczego zróżnicowano je pod względem jakościowym. Szerszego wyjaśnienia wymagałyby także zmienne opisujące natężenie połowów. Jego brak uniemożliwia właściwą interpretację wyników na ich podstawie.

Wczytanie danych

Dane źródłowe dla naszej analizy zostały umieszczone w pliku sledzie.csv. Poniższy blok kodu prezentuje wczytanie danych z pliku źródłowego i zapisanie ich jako data frame sledzie.

sledzie <- read.csv("~/R/workspace/Rcwicze/sledzie.csv",  na.strings="?")

Poniżej dziesięć pięrwszych wierszy z zestawu danych. Zauważmy, że w pierwszym wierszu tabeli występuje brakująca wartość atrybutu chel2.

knitr::kable(head(sledzie,n = 5))
X length cfin1 cfin2 chel1 chel2 lcop1 lcop2 fbar recr cumf totaln sst sal xmonth nao
0 23.0 0.02778 0.27785 2.46875 NA 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
1 22.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
2 25.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
3 25.5 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8
4 24.0 0.02778 0.27785 2.46875 21.43548 2.54787 26.35881 0.356 482831 0.3059879 267380.8 14.30693 35.51234 7 2.8 
str(sledzie)
## 'data.frame':    52582 obs. of  16 variables:
##  $ X     : int  0 1 2 3 4 5 6 7 8 9 ...
##  $ length: num  23 22.5 25 25.5 24 22 24 23.5 22.5 22.5 ...
##  $ cfin1 : num  0.0278 0.0278 0.0278 0.0278 0.0278 ...
##  $ cfin2 : num  0.278 0.278 0.278 0.278 0.278 ...
##  $ chel1 : num  2.47 2.47 2.47 2.47 2.47 ...
##  $ chel2 : num  NA 21.4 21.4 21.4 21.4 ...
##  $ lcop1 : num  2.55 2.55 2.55 2.55 2.55 ...
##  $ lcop2 : num  26.4 26.4 26.4 26.4 26.4 ...
##  $ fbar  : num  0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 0.356 ...
##  $ recr  : int  482831 482831 482831 482831 482831 482831 482831 482831 482831 482831 ...
##  $ cumf  : num  0.306 0.306 0.306 0.306 0.306 ...
##  $ totaln: num  267381 267381 267381 267381 267381 ...
##  $ sst   : num  14.3 14.3 14.3 14.3 14.3 ...
##  $ sal   : num  35.5 35.5 35.5 35.5 35.5 ...
##  $ xmonth: int  7 7 7 7 7 7 7 7 7 7 ...
##  $ nao   : num  2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 2.8 ...

Wstępne przetwarzanie danych

Poniższy blok kodu generuje statystyki opisowe źródłowego zbioru danych. W obliczeniach pomijamy zmienne:

  • X, która jako zmienna porządkowa nie wnosi nic do naszej analizy,
  • xmonth, ograniczymy się do sprawdzenia, czy zawiera jedynie poprawne wartości (Liczby całkowite 1,2,..,12).

Dla pozostałych zmiennych odnotujmy następujące obserwacje:

  • Atrybut fbar, oraz cumf zawierają się w przedziale (0,1).
  • nao jest jedyną zmienną przyjmującą wartości ujemne.
  • Różnica między najwyższą temperaturą przy powierzchni wody (sst) na przestrzeni 60 lat wynosiła 2°C.
  • Najmniejszą zmiennością cechuje się zasolenie wody. Cały jej zakres wynosi (0,2 ppt)
  • Wartości puste występują dla siedmiu zmiennych. Około 1600 przypadków dla każdej z nich.
summary(sledzie[-c(1,15)])
##      length         cfin1             cfin2             chel1       
##  Min.   :19.0   Min.   : 0.0000   Min.   : 0.0000   Min.   : 0.000  
##  1st Qu.:24.0   1st Qu.: 0.0000   1st Qu.: 0.2778   1st Qu.: 2.469  
##  Median :25.5   Median : 0.1111   Median : 0.7012   Median : 5.750  
##  Mean   :25.3   Mean   : 0.4458   Mean   : 2.0248   Mean   :10.006  
##  3rd Qu.:26.5   3rd Qu.: 0.3333   3rd Qu.: 1.7936   3rd Qu.:11.500  
##  Max.   :32.5   Max.   :37.6667   Max.   :19.3958   Max.   :75.000  
##                 NA's   :1581      NA's   :1536      NA's   :1555    
##      chel2            lcop1              lcop2             fbar       
##  Min.   : 5.238   Min.   :  0.3074   Min.   : 7.849   Min.   :0.0680  
##  1st Qu.:13.427   1st Qu.:  2.5479   1st Qu.:17.808   1st Qu.:0.2270  
##  Median :21.673   Median :  7.0000   Median :24.859   Median :0.3320  
##  Mean   :21.221   Mean   : 12.8108   Mean   :28.419   Mean   :0.3304  
##  3rd Qu.:27.193   3rd Qu.: 21.2315   3rd Qu.:37.232   3rd Qu.:0.4560  
##  Max.   :57.706   Max.   :115.5833   Max.   :68.736   Max.   :0.8490  
##  NA's   :1556     NA's   :1653       NA's   :1591                     
##       recr              cumf             totaln             sst       
##  Min.   : 140515   Min.   :0.06833   Min.   : 144137   Min.   :12.77  
##  1st Qu.: 360061   1st Qu.:0.14809   1st Qu.: 306068   1st Qu.:13.60  
##  Median : 421391   Median :0.23191   Median : 539558   Median :13.86  
##  Mean   : 520367   Mean   :0.22981   Mean   : 514973   Mean   :13.87  
##  3rd Qu.: 724151   3rd Qu.:0.29803   3rd Qu.: 730351   3rd Qu.:14.16  
##  Max.   :1565890   Max.   :0.39801   Max.   :1015595   Max.   :14.73  
##                                                        NA's   :1584   
##       sal             nao          
##  Min.   :35.40   Min.   :-4.89000  
##  1st Qu.:35.51   1st Qu.:-1.89000  
##  Median :35.51   Median : 0.20000  
##  Mean   :35.51   Mean   :-0.09236  
##  3rd Qu.:35.52   3rd Qu.: 1.63000  
##  Max.   :35.61   Max.   : 5.08000  
##

Przetwarzenie wartości pustych

Rzut okiem na pierwsze 5 wierszy wczytanego pliku sugeruje, że wartości puste nie dominują w danych surowych. Sprawdźmy jaki wpływ na rozmiar danych będzie miało ich usunięcie.

tb_sledzie <- tbl_df(sledzie) %>%
    filter(!is.na(chel1), !is.na(chel2), !is.na(cfin1) , !is.na(cfin2) , !is.na(lcop1) , !is.na(lcop2), !is.na(sst))

Liczba obserwacji data frame’u sledzie po usunięciu wartości pustych to 42488. Stanowi to 81% rozmiaru początkowego.

Analiza atrybutów

Poniższy kod wylicza ile unikalnych wartości posiada każdy z atrybutów.

unique_tb_sledzie <- tb_sledzie %>% 
  summarise(length = n_distinct(length), cfin1 = n_distinct(cfin1), cfin2 = n_distinct(cfin2), chel1 = n_distinct(chel1), chel2 = n_distinct(chel2), lcop1 = n_distinct(lcop1),   lcop2 =n_distinct(lcop2), fbar = n_distinct(fbar), recr = n_distinct(recr), cumf = n_distinct(cumf), totaln = n_distinct(totaln), 
            sst = n_distinct(sst), sal = n_distinct(sal), nao = n_distinct(nao)) 

unique_tb_sledzie_t <- as.data.frame(cbind( as.matrix(names(unique_tb_sledzie)),  t(unique_tb_sledzie)))



ggplot(unique_tb_sledzie_t, aes(x=reorder(V1, -as.integer(V2)), y = V2)) + 
    geom_bar(stat = "identity") + 
    labs(title ="liczba unikalnych wartości każdego z atrybutów", x = "nazwa zmiennej", y="liczba unikalnych wartości")

Największą różnorodnością cechuje się zmienna length. Chociaż 55 unikalnych wartość przy blisko pięćdziesięciu tysiącach obserwacji wydaje się liczbą znikomą dla atrybutu ciągłego to zauważmy, że minimalna długość sledzia w naszym zbiorze danych wynosi 19cm natomiast maksymalna 32.5cm. Przy założeniu, że pomiarów dokonywano z dokładnością do 0.5 centymetra, zmienna length powinna mieć 27 unikalnych wartości, natomiast przy dokładności 0,01cm- 135.

Długość (length)

spójrzmy na częstości zmiennej length

table(tb_sledzie$length)
## 
##   19 19.5   20 20.5   21 21.5   22 22.2 22.4 22.5   23 23.2 23.3 23.4 23.5 
##    2    5   23   69  173  385  791    1    1 1254 2172    1    6    1 2712 
## 23.6 23.7 23.8 23.9   24 24.1 24.2 24.4 24.5 24.7 24.8 24.9   25 25.1 25.3 
##    1    2    6    1 3603    1    5    2 4345    6    4    3 4763    2    1 
## 25.4 25.5 25.7 25.8   26 26.1 26.3 26.5 26.8 26.9   27 27.3 27.5   28 28.5 
##    3 4868    2    1 4672    1    1 4192    1    1 3265    2 2381 1498  782 
## 28.6   29 29.5 29.8   30 30.5   31 31.5   32 32.5 
##    1  301  113    1   24   21    9    3    2    2

Zwróćmy uwagę, że tylko dla pojedyńczych obserwacji wynik pomiaru zapisano z dokładnością do 1mm. Przyjmijmy zatem, że pomiaru dokonano z dokładnością do 5mm. Dla uspójnienia wyników zaokrąglimy części dziesiętne z dokładnością do 0,5cm.

tb_sledzie <- tb_sledzie %>%
  mutate(length = round(length*2)/2)

table(tb_sledzie$length)
## 
##   19 19.5   20 20.5   21 21.5   22 22.5   23 23.5   24 24.5   25 25.5   26 
##    2    5   23   69  173  385  792 1255 2173 2722 3616 4353 4772 4874 4674 
## 26.5   27 27.5   28 28.5   29 29.5   30 30.5   31 31.5   32 32.5 
## 4193 3267 2383 1498  783  301  113   25   21    9    3    2    2

Histogram tak uproszczonej zmiennej length prezentuje wykres poniżej. Czerwoną linią zaznaczono krzywą funkcji gęstości rozkładu normalnego o średniej i wariancji takich jak zmienna length.

sledzie_hist <- ggplot(tb_sledzie, aes(x = length)) + geom_bar() + theme_light()

n <- nrow(tb_sledzie)
mean <- mean(tb_sledzie$length)
sd <- sd(tb_sledzie$length)
binwidth = 0.5 # passed to geom_histogram and stat_function
set.seed(1)
df <- data.frame(x = rnorm(n, mean, sd))

sledzie_hist + stat_function(fun = function(x) dnorm(x, mean = mean, sd = sd) * n * binwidth , color = "darkred", size = 1)

Wykres wygląda na nieco spłaszczony w stosunku do krzywej. Lewy ogon naszego histogramu jest wyraźnie grubszy niż prawy. Podobne wnioski można wyciągnąć z wzrokowej oceny wykresu kwantyl-kwantyl dla zmiennej length (poniżej). Nasycenie lewego ogona może być spowodowane przez zjawisko karłowacenia. Na podstawie danych nie jesteśmy jednak w stanie określić jak kształtował się rozmiar śledzi w kolejnych latach.

ggplot(tb_sledzie, aes(sample = length)) + stat_qq() + stat_qq_line() + theme_light()

Dane na temat dostępności plankotnu

Poniższy blok kodu generuje wykres pudełkowy dla zmiennych opisujących zagęszczenie planktonu w miejscach gdzie dokonywano połowów.

boxplot(tb_sledzie[,3:8], xlab ="gatunek planktonu", ylab="dostępność planktonu", main = "Rozkłady zagęszczenia rozważanych gatunków planktonu")

Opis zbioru przeznaczonego do analizy nie zawiera informacji na temat metodologii pomiaru zagęszczenia. Jednak wykresy pudełkowe dla Calanus finmarchicus (cfin1 i cfin2) pokazuję iż zdecydowana większość pomiarów jest skupiona wokół zera. W szczególności cfin1, w porównaniu z innymi gatunkami widłonogów, niemal nie występował na łowiskach gdzie dokonano obserwacji.

hist_cfin1 <- ggplot(tb_sledzie, aes(cfin1)) + geom_histogram(bins = 38) + theme_light() + labs(title="Calanus finmarchicus gat. 1", x="", y = "")
hist_cfin2 <- ggplot(tb_sledzie, aes(cfin2)) + geom_histogram(bins = 20) + theme_light() + labs(title="Calanus finmarchicus gat. 2", x="", y = "")
grid.arrange(hist_cfin1, hist_cfin2, nrow=1, top = "Histogram zagęszczenia", left ="Liczba obserwacji")

Lepiej wygląda zagęszczenie Calanus helgolandicus. Obserwowane wartości są wyraźnie większe od zera. Choć na wykresie pudełkowym duże wartości (60, 80) są oznaczone jako odstające, histogramy pokazują, że wystąpiły blisko lub ponad tysiąc krotnie.

hist_chel1 <- ggplot(tb_sledzie, aes(chel1)) + geom_histogram(bins = 80) + theme_light() + labs(title="Calanus helgolandicus gat. 1", x="", y = "")
hist_chel2 <- ggplot(tb_sledzie, aes(chel2)) + geom_histogram(bins = 60) + theme_light() + labs(title="Calanus helgolandicus gat. 2", x="", y = "")
grid.arrange(hist_chel1, hist_chel2, nrow=1, top = "Histogram zagęszczenia", left ="Liczba obserwacji")

Histogramy dla pomiarów zagęszczenia widłonogów również pokazały, iż wartości przedstawione na wykresie pudełkowym jako odstające występowały stosunkowo często w stosunku do innych wyników.

hist_lcop1 <- ggplot(tb_sledzie, aes(lcop1)) + geom_histogram(bins = 80) + theme_light() + labs(title="Widłonogi gat. 1", x="", y = "")
hist_lcop2 <- ggplot(tb_sledzie, aes(lcop2)) + geom_histogram(bins = 80) + theme_light() + labs(title="Widłonogi gat. 2", x="", y = "")
grid.arrange(hist_lcop1, hist_lcop2, nrow=1, top = "Histogram zagęszczenia", left ="Liczba obserwacji")

Poniżej macierz korelacji zmiennych opisujących zagęszczenie planktonu i długości śledzia. Wyznaczone współczynniki pokazują, że długość śledzia nie jest liniowo skorelowana z zagęszczeniem planktonu. Liniowa korelacja występuje natomiast pomiędzy zmiennymi lcop1 oraz chel1. Wypływa stąd wniosek, że calanus helgolandicus jest głównym gatunkiem widłonogów występujących na obszarze połowów objętych badaniem i jego występowanie decyduje o zagęszczeniu widłonogów w ogóle.

widlonogi_cor_mat <- cor(tb_sledzie[,2:8])

corrplot(widlonogi_cor_mat, method="color", addCoef.col = "gray")

Powyższe wnioski potwierdzają także obserwacje diagramów punktowych dla omówionych par zmiennych.

chel1_lcop1_point <- ggplot(tb_sledzie, aes(y=chel1, x = lcop1)) + geom_point() + theme_light() + geom_smooth(color="red", method = "lm") + labs(x = "widłonogi 1. gat.", y = "calanus helgolandicus 1. gat.")
chel2_lcop2_point <- ggplot(tb_sledzie, aes(y=chel2, x = lcop2)) + geom_point() + theme_light() + geom_smooth(color="red", method = "lm") + labs(x = "widłonogi 2. gat.", y = "calanus helgolandicus 2. gat.")

cfin1_lcop1_point <- ggplot(tb_sledzie, aes(y=cfin1, x = lcop1)) + geom_point() + theme_light() + labs(x = "widłonogi 1. gat.", y = "calanus finmarchicus 1. gat.")
cfin2_lcop2_point <- ggplot(tb_sledzie, aes(y=cfin2, x = lcop2)) + geom_point() + theme_light() + labs(x = "widłonogi 2. gat.", y = "calanus finmarchicus 2. gat.")

grid.arrange(chel1_lcop1_point, chel2_lcop2_point, nrow=1, top = "Związek pomiędzy zagęszczeniem calanus helgolandicus a zagęszczeniem wszystkich widłonogów")

grid.arrange(cfin1_lcop1_point, cfin2_lcop2_point, nrow=1, top = "Związek pomiędzy zagęszczeniem calanus finmarchicus a zagęszczeniem wszystkich widłonogów")

length_lcop1_point <- ggplot(tb_sledzie, aes(y=length, x = lcop1)) + geom_point() + geom_smooth(color = "red", method = "lm") + theme_light() + labs(x = "widłonogi gat. 1.", y = "lenght", title="Diagram punktowy: korelacja między długością śledzia, a zagęszczeniem widłonogów gat 1.")

length_lcop1_point

Na podstawie powyższych obserwacji w dalszej analizie pominiemy zmienne opisujące zagęszczenie planktonu za wyjątkiem lcop1, który jest najsilniej skorelowany liniowo z długością śledzia.

Badanie korelacji liniowej.

W tej części zaprezentujemy analizę korelacji między zmienną length a pozostałymi atrybutami. Pochylimy się także nad zależnościami pomiędzy pozotałymi atrubutami opisowymi.

Macierz korelacji dla tabeli tb_sledzie, po usunięciu zmiennej porządkowej i zredukowaniu zmiennych opisujących zagęszczenie planktonu do lcop1, prezentuje poniższy wykres. Łatwo zauważyć, że najsilniej skorelowana z atrybutem length jest zmienna sst opisująca temperaturę wody. Nie jest to jednak wartość, która sugerowałaby zależność liniową.

Współczynnik korelcji pearsona wydaje się też odkrywać pewien wpływ na długość zmiennych fbar oraz nao.

all_cor_mat <- cor(tb_sledzie[,-c(1,3,4,5,6,8)])

corrplot(all_cor_mat, method="color", addCoef.col = "gray")

Interesujący jest brak jakiejkolwiek korelacji między zmienną xmonth, a pozostałymi atrybutami.

Silna zależność liniowa występuje pomiędzy zmiennymi fbar i cumf opisującymi odpowiednio natężenie połowów w regionie oraz łaczne roczne natężenie połowów w regioniu. Może to oznaczać, że natężenie połowów na danym obszarze nie zmienia się istotnie w ciągu roku. Brak zjawiska sezonowości w natężeniu połowów tłumaczy także brak korelacji między zmienną xmonth a pozostałymi atrubutami. Nieco mniej ale nadal dość istotnie powiązane są zmienne fbar i cumf ze zmienną totaln. Związek między łacznym rocznym natężeniem połowów, a liczbą śledzi złowionych w ramach danego połowu jest silniejszy niż między natężeniem połowów w regionie. Wydaje się to dziwne i nie potrafimy podać sensownej interpretacji tego zjawiska.

Podsumowanie analizy korelacji

Wykres poniżej prezentuje macierz korelacji dla zmiennych, które cechują się najwyższym bezwzglednym współczynnikiem korelacji pearsona z długością śledzia. Zmienne te nie są jednak niezależne między sobą. Dominujący wpływ na pozostałe czynniki wydaje się mieć temperatura wody (sst). W szczególności atrybut nao, czyli oscylacja północno atlantycka, którą pominiemy w dalszej analizie.

red_all_cor_mat <- cor(tb_sledzie[,c(2,7,9,13,16)])

corrplot(red_all_cor_mat, method="color", addCoef.col = "gray")

Zależność między zmianami przeciętnej długości śledzia a temperaturą.

Dane nie są opatrzone datami. Zakładamy jednak, że kolejne obserwacje zapisane w zbiorze były rejestrowane chronologicznie. Na diagramie punktowym zaznaczyliśmy linię trędu pokazującą przeciętną długość złowionych śledzi. Nasycenie barwą niebieską prezentuje temperaturę zarejestrowaną dla danego połowu. Błękit oznacza temperaturę najwyższą, ciemny granat najniższą. Na wykresie widać że czerwona linia pnie się do góry na tle zgrupowania ciemnych kropek, natomiast opada wraz z przewagą jasnych kropek w późniejszym okresie. Potwierdza to istotny wpływ temperatury na zmianę przeciętnej długości śledzia.

p_mean_len_temp_time <- ggplot(tb_sledzie, aes(X, length, color = sst)) + geom_point() + geom_smooth(method = "auto", color = "red") + theme_light() + labs(title="zmiana długości śledzi")


ggplotly(p_mean_len_temp_time)

Predykcja

W niniejszym rozdziale zaprezentujemy metody uczenia maszynowego: podział zbioru na uczący i testowy, stratyfikację danych oraz miary oceny regresji. Wykorzystamy do tego bibliotekę caret.

Podział zbioru danych na zbiór uczący i testowy.

Poniższy blok kodu rozpoczyna usunięcie wartości odstających, które ze względu na zbyt małą liczebność uniemożliwiają skuteczną stratyfikację zbioru danych. Wartości odstające wykorzystamy później dla sprawdzenia jak poradzi sobie z nimi wytrenowany przez nas model regresji.

to_train_tb_sledzie <- tb_sledzie %>%
    select(length, lcop1, fbar, sst) %>%
    filter(length > 19 & length < 31.5)

sledzie_outliers <- tb_sledzie %>%
    select(length, lcop1, fbar, sst) %>%
    filter(length <= 19 | length >= 31.5)

set.seed(11)
inTraining <- createDataPartition(
                      y = to_train_tb_sledzie$length,
                      p = 0.7,
                      list = FALSE)

training <- to_train_tb_sledzie[inTraining,]
testing <- to_train_tb_sledzie[-inTraining,]
ctrl <- trainControl(
    # powtórzona ocena krzyżowa
    method = "repeatedcv",
    # liczba podziałów
    number = 5,
    # liczba powtórzeń
    repeats = 5)

Wybór i dopasowanie modelu, algorytm Random Forest

set.seed(11)
fit_1 <- train(length ~ .,
             data = training,
             method = "rf",
             importance=T,
             trControl = ctrl,
             ntree = 10)
## note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
set.seed(11)
fit_ntree_30 <- train(length ~ .,
             data = training,
             importance=T,
             method = "rf",
             trControl = ctrl,
             ntree = 30)
## note: only 2 unique complexity parameters in default grid. Truncating the grid to 2 .
rfGrid <- expand.grid(mtry = 1:3)
gridCtrl <- trainControl(
    method = "repeatedcv",
    number = 5,
    repeats = 5)

set.seed(11)
fit_mtry_tune <- train(length ~ .,
             data = training,
             method = "rf",
             metric = "RMSE",
             importance=T,
             trControl = gridCtrl,
             tuneGrid = rfGrid,
             ntree = 5)
set.seed(11)
fit_tune <- train(length ~ .,
             data = training,
             method = "rf",
             metric = "RMSE",
             importance=T,
             trControl = gridCtrl,
             tuneGrid = expand.grid(mtry = 3),
             ntree = 10)

Poniższe wykresy prezentują kolejne kroki dopasowywania modelu mierzonego przy pomocy pierwiastka błędu średniokwadratowego dla różnych wartości parametru mtry. Pierwszy model został zaproponowany automatycznie. kolejny pokazuje nieznaczą poprawę dopasowania po zwiększeniu parametru ntree, natomiast ostatnio pokazuje badanie jakości dopasowanie dla całkowitych wartości parametru mtry z przedziału [2,20].

ggp1 <- ggplot(fit_1) + theme_light() + labs(title = "ntree = 10", y="", x="mtry")
ggp2 <- ggplot(fit_ntree_30) + theme_light() + labs(title = "ntree = 30", y="", x="mtry")
ggp3 <- ggplot(fit_mtry_tune) + theme_light() + labs(title = "ntree = 5", y="", x="mtry")

grid.arrange(ggp1, ggp2, ggp3, nrow=1, left="RMSE (Repeated Cross-Validation)", top = "Dopasowanie modelu")

kable(fit_1$results[,1:4])
mtry RMSE Rsquared MAE
2 1.190307 0.4792417 0.9431709
3 1.190645 0.4789517 0.9432560 
kable(fit_mtry_tune$results[,1:4])
mtry RMSE Rsquared MAE
1 1.190543 0.4790450 0.9434371
2 1.190518 0.4790859 0.9433590
3 1.190959 0.4786796 0.9435271 
kable(fit_ntree_30$results[,1:4])
mtry RMSE Rsquared MAE
2 1.190086 0.4794358 0.9429409
3 1.190502 0.4790734 0.9431261

Ocena ważności parametru

plot(varImp(fit_tune))

Analiza ważności atrybutów pokazuje, że najistotniejsza dla opisu zmienności długości śledzia jest temperatura przy powierzchni wody. Pozostałe parametry zostały ocenione jako zdecydowanie mniej istotne.

Na koniec pokażemy jak wytrenowany przez nas model Random Forest poradził sobie z danymi testowymi.

rfPredicts <- predict(fit_tune, newdata = testing)

rf_outliresPredicts <- predict(fit_tune, newdata=sledzie_outliers)

RMSE(rfPredicts, testing$length)
## [1] 1.175606